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Abstract. We report on modelling in stellar astrophysics with the ANTARES code. First, we 
describe properties of turbulence in solar granulation as seen in high-resolution calculations. 
Then, we turn to the first 2D model of pulsation-convection interaction in a cepheid. We discuss 
properties of the outer and the Hell ionization zone. Thirdly, we report on our work regarding 



CO models of semiconvection in the context of stellar physics. 
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!■ In this paper we present applications of the ANTARES code to various problems in 
stellar physics. ANTARES is a radiation-hydrodynamics code. Its core functionality has 



I 

5_i ■ recently been described by Muthsam et al. (2010) Using this core functionality, we first 



discuss properties of turbulence in solar granulation, studied in extremely high resolution 
d '. in 3D. 

We have in the meanwhile extended the core functionality of ANTARES, and the ap- 
-H ■ plications discussed in the next two chapters make use of these extensions. So, we have 
included a spherical and optionally radially moving grid for the study of the interaction 
of radial stellar pulsation with convection. We present a first description of the pulsation- 
convection interaction of a 2D model of a cepheid with realistic equation of state and grey 



^v^j . radiative transfer. For the importance of studying the pulsation-convection interactions 



consult, for example, the reviews |Buchler (1997)] and |Brichler (2009) | Consider that since 



. the pioneering work due to Christy, Cox and Kippenhahn in the 1960's improvement in 
nonlinear stellar pulsation modelling has basically been due to better microphysics and to 
t-H . some extent to improvement in numerics. Convection has, however, (practically) always 
been included by a mixing-length type approach where already the basic form of the 
equations may be prone to some doubt and where additionally a number of parameters 
needs to be set in a way which can barely considered to be really convincing. Problems 



probably stemming from such shortcomings are discussed, e.g., in Buchler (1997)[ For a 
recent discussion of a specific uncertainty which arises from traditional convection mod- 
elling in pulsating stars (namely, excitation of double mode pulsations or lack thereof) 
see 



Smolec fc Moskalik (2009) 



Subsequently, we turn our attention towards modelling of semiconvection in the con- 
text of stellar physics. Although it is textbook knowledge that a faithful description of 
semiconvection is highly important for properly dealing with the late stages of stellar evo- 
lution, but little work has been undertaken with regard to multidimensional modelling 
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in the stellar context. The only major numerical studies are due to Merryfield (1995) 



and Biello (2001) See also the results due to Bascoul (2007) We present here our first 
2D results discussing, in particular, the properties of diffusive interfaces. These results 
mainly refer to Boussinesq models, although we report also on compressible models and 
discuss their connection with the Boussinesq case. For enabling a more extensive investi- 
gation of compressible models it is necessary to implement a method which allows large 
time-steps, restricted by the macroscopic rather than the sound velocity. In this connec- 
tion, an implementation of the method of Kwatra et al.(2009) is in an advanced stage of 
development. 



2. High-resolution solar granulation simulations 

In order to investigate the turbulent state of solar granulation we have performed hy- 
drodynamic simulations with extremely high resolution. While they are basically similar 
to the ones described in Muthsam et al. (2010)| they differ in that the highest resolution 
is achieved by making use of two grid refinement zones (stacked within each other) in- 




Figure 1. Solar granulation at two instances of time. Horizontal extent of the domain shown: 
1 Mm. The grey surface is an isosurface for T = 6000 K. The vortex tubes are made visible by 
plotting isosurface for a (large) value of the modulus of the gradient of a suitably normalized 
pressure. A comparison of the upper and lower figure shows the effect of increased degree of 
turbulence when applying higher resolution. 
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stead of one, and the resolution in the smaller one is considerably finer than in the old 
run. In addition, the old run referred to an exploding granule, whereas we concentrate 
here, in the highest resolution, on a strong downdraft surrounded by normal granules in 
order to figure out to what extent our previous results apply here as well. 

More specifically, the basic computational domain spans now 3.68 x 6 x 6 Mm 3 (the 
first coordinate is the vertical) with a cell size of 10.8 x 22.2 x 22.2 km 3 . The first grid 
refinement zone has an extension of 1.9 x 2.5 x 2.8 Mm 3 and a grid spacing of 5.4 x 7.4 x 
7.4 km 3 . The second and therefore finest grid refinement zone contains the downdraft 
mentioned above and parts of the surrounding granules. It has an extent of 1 x 1.2 x 
1.2 Mm 3 with a cell size of 2.7 x 3.7 x 3.7 km 3 . 

Using the high-resolution calculations, we want to discuss the role of resolution in the 
representation of the turbulent state of the solar granulation layer. In earlier simulations 
of an exploding granule Muthsam et al. (2010) a host of vortex tubes of small diameter 
has shown up, located preferentially near the granular lanes and the downdrafts. The 
questions popping up are now whether a similar large number of vortex tube shows 
up even in more normal granulation and whether these previous calculations were fully 
resolved; they used a grid spacing of 7.1 x 9.8 x 9.8 km 3 , i.e. worse than even our present 
first refinement. 

In order to answer the question at what resolution the turbulent field of solar granu- 
lation is (possibly) fully resolved, we compare in figure Q] two snapshots of our highest 
resolution subdomain. The first figure represents the state shortly after we have turned 
on the second grid refinement and thus still have the effective resolution due to the first 
refinement, whereas the next snapshot is taken 153 seconds later. Due to increased res- 
olution, the number and intensity of the vortex tubes has drastically increased. A closer 
study of the properties of this special sort of turbulence encountered mainly in granular 
downdrafts is under way. 



3. Cepheid pulsation and convection: a 2D model 

Initial condition and computational domain. We model a cepheid with an effective 
temperature of 5125K, luminosity L = 912.797L Q , mass M — 5.0Af o , hydrogen content 
X = 0.7 and metallicity Z = 0.01. Our computational domain reaches from 4000 K to 
320.000 K, thus the outer 42% of the star are computed. 

This domain is equipped with a polar grid. In radial direction N x — 510 grid points 
plus 4 ghost cells at each boundary are used. The r-range covers r £ [r top , r{, ot } , where r\, ot 
is fixed at 15.5 Gm and r top varies with time from about 26.2 Gm to 27.4 Gm. The grid is 
stretched in radial direction by a factor q. The mesh sizes are Ari + \=qAri varying from 
A7'o = 0.046 Mm at the top to Ar^^-i = 12 Mm at the bottom. In angular direction 
there are 800 grid points, the distance riAip between two adjacent gridpoints in the 
H-ionisation zone is 5.8 Mm and 5.6 Mm in the He-ionisation zone. The corresponding 
aspect ratios are 1 : 2.5 and 1 : 0.6. 

Fluxes and pulsation in the Hell-ionisation zone. The work integral can be decomposed 
into an average part PdVo — d r uo(jp/ p) and a perturbational part defined as PdV per t = 
d r u r (p/p). Here, uo = I r /p and u' r = u r — uq. I = (I r , J v ) is the momentum vector 
and u = (u r , u v ) the velocity. I r denotes the horizontal average of the radial momentum 
component. 

While PdVp er t varies locally from — 6.e8 to +5.e8 in the Hell convection zone and 
is greater than PdVo £ [— 8.e7, +5.e7], the horizontal average is considerably smaller, 
PdVp er t € [-9.e6, +7.e6], compared to PdV ~ PdV . 

The averaged convective flux varies from —3% to +8% of the input energy flux at the 
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Figure 2. Convective flux in the course of contraction (Hell convection zone) 




Figure 3. Horizontal averages of convective flux (upper figure) and work integral (lower figure) 

lower boundary, the mean value during one period being +4%. Locally values of up to 
1.4ell corresponding to 30% can be reached. The averaged kinetic flux varies from —3% 
to +0.2% the mean value is —2%. The data are from 12 consecutive periods and the 
same flux pattern can be observed in each though the extent may vary. In the expanded 
state one observes two convection centers: a new one forming at the top and an older 
one which is the remainder of the downwards moving and slowly disapearing plume of 
the previous period. At contraction the two centers are very close to each other. These 
centers of convective activity lead also to the two (occasionally more) stripes situated 
atop of each other and visible, for example, in the quantities depicted in figure [3] The 
centers of activity can also directly be seen in figure [5] 

Remarks on the H-ionisation zone. The upper boundary of our computational do- 
main is situated at an optical deph of a few thousandths. Figure 2] shows a part (in 
4>— direction) of the upper zone. Color coding is for momentum densities (bright resp. 
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Figure 4. Radial momentum in H-ionisation zone. The bright line denotes the optical depth 

t = 100 

red is downwards). The convective motions are very vigorous. We have observed Mach 
numbers for inflow of up to 4 and for outflow of up to 3. Remarkably, the lower boundary 
of the downdrafts often coincides the line of maximum temperature gradient (just above 
the bright line in figure|4j These results could only be achieved with grid refinement. The 
refined grid spacing is one third of the coarse grid spacing in the radial and one fourth 
in the angular direction. 



4. Numerical simulation of semiconvection 

When the convective transport does not depend on the temperature gradient only 
but also on an additional scalar field like salt (in the ocean) or helium (in stars), then 
double-diffusive convection can occur. This name implies that two diffusion coefficients 
(kt and ks/khc) influence the fluids motion. In the oceans double-diffusive convection 
can form salt-fingers, while the well-known Latte Macchiato layers are a result of another 
double-diffusive phenomenon, which in astrophysics is called semi-convection. 

In the evolution of stars semi-convection plays an important role. When cold hydrogen- 
rich matter is stratified above hot helium-rich matter, double-diffusive layering can occur 
under certain circumstances (Ledoux stable, V M > V — V oc j)- The main question, which 
arises is the influence of the double-diffusive mixing processes on the evolution of high 
mass stars (M > 15Mq) around main sequence turnoff. 

A stable layered structure in semi-convection zones (cf. Latte Macchiato), separated by 
diffusive boundary layers, is assumed (|Huppert fc Moore (1976)||Spruit(T992) ). In partic- 



ular, the mass and heat transport through these diffusive boundaries are of major interest 
for explaining mixing time scales and the life span of the entire semi- convective region. 
There are several unknowns when describing single layers. Especially, the thickness of 
a single layer as function of thermal and saline Rayleigh number (Rax, Ras) an d the 
superadiabatic gradient is not well understood. A comprehensive hydrodynamical model, 



based on Spruit(1992) and Muthsam et al. (1999) (see also Muthsam et al.(1995)) was 



derived for the double diffusive convection case to tackle these uncertainties. 



In Zaussinger & Spruit(2010) numerical simulations in 2D of semi-convection, based 
on compressible and incompressible formulations, have been performed for a wide initial 
parameter range of the Prandtl number a = v/kt, the Lewis number r = ks/kt, the 
modified Rayleigh number Ra„ — Rax (J and the stability parameter R p = Ras / Rax- 
It is proven there that simulations done in the Boussinesq approximation are directly 
comparable with fully compressible fluids as long as the layer heights are small enough 
(Hp > > H) . The aim of this study was the verification and extent of validity of existing 
power laws in the form Nut = aRa T , where Nut is the thermal Nusselt number. An 
extrapolation to the stellar parameter space was done subsequently. 



In Zaussinger(2010) the parameter dependence of the fluxes is analysed. Single and 
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double layer simulations were at the focus of that study. While the diffusion associated 
dimensionless numbers (a, r) can be compared directly between compressible and in- 
compressible fluids, other physical quantities like the fluxes (Ft, Fsu e ) or the Rayleigh 
numbers (Ra,T, Ras) need a special treatment. The thermal Nusselt number, which com- 
pares the diffusive heat flux to the total heat flux, has to be corrected by the adiabatic 
stratification flux F ac i, which is by definition not present in incompressible fluids mod- 
elled by the Boussinesq approximation. Neglecting this difference would lead to wrong 
comparisons in simulations with flat temperature gradients. 
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Figure 5. Nusselt numbers as function of time. The formation of a semi-convective zone (3) is 
the consequence of a diffusive phase (1) and a oscillatory phase (2). For the specified parameter 
space the convective rolls are established within 12 turnover cycles. 

Linear stability analysis shows that semi-convection in an oscillatory instability. This 
could be observed in all our simulations. Figure [5] illustrates the development of a semi- 
convection cell in terms of the Nusselt number when starting from a linear stratification. 
Nusselt numbers in the range 1 < Nut < 5 indicate very diffusive regimes, where the 
total heat flux is dominated by conduction. High values correspond to more convective 
regimes. A diffusive phase (1) is followed by an oscillatory phase (2). The fully evolved 



semi-convective role (3) leads to the expected (Spruit(1992) ) step like structure in helium 
and temperature at the upper and lower boundaries. With a suitable choice of parame- 
ters, a preassigned step in the middle of the domain is found to be stable for many eddy 
turnover times. 

While adoption of linear stratification as initial condition leads to very long simulation 
times, a 'direct initial step stratification' is considered. It has been shown that both 
approaches lead to the same results. Strong damping mechanisms in the parameter space 
R p > 2 and Le > fO -2 lead to very long simulation times. By starting from a step like 
stratification these could be reduced significantly and a formerly inaccessible parameter 
range came into reach. 

The verification of theoretical and experimental power laws under consideration for the 
stellar parameter space was the main focus of this study. It turned out that several ex- 



isting 'fitting formulas' ( |Castaing et al. (1989)[ |Spruit(1992")| |Niemela et al. (2000) [ ) de 



scribe an upper limit. Even the flux relation, 

Nu s = r 1/2 Nu T (4.1) 
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could be verified for semi-convective stable layers (see figure . Equipped with these 
results an extension to the stellar parameter regime was considered. The main problem of 
existing ID stellar evolution codes is the lack of information about to the semi-convective 
mixing efficiency. The superadiabatic gradient V — V a d is not known in advance and has 
to be estimated. The semi-convective zone in a stellar model is limited in extent and 
consequently the evolution of the star on small time-scales does not depend much on the 
way it is calculated. But the additional diffusion in helium could be important for later 
evolutionary stages. This leads to the assumption that the heat flux could be considered 



as known rather than the real temperature gradient V. Based on Spruit(1992) an implicit 



function was derived, relating the temperature gradient, the Rayleigh number, the fluxes 
and the the vertical extent of the semi-convective zone to each other. 




R =1.15 Pr=1.0 Le=0.10 
R^=1.50 Pr=1.0 Le=0.10 
FT =2.00 Pr=1.0 Le=0.10 
Ro=5.00 Pr=1.0 Le=0.10 
Nu T = sqrt(Le=0.1)-Nu s 
R =1.15 Pr=0.1 Le=0.1 - 
R„=1.50Pr=0.1 Le=0.1 
FT =2.00 Pr=0.1 Le=0.1 - 
rr=5.00 Pr=0.1 Le=0.1 
Nu T = sqrt(Le=0.01)'Nu s 
R =1.15 Pr=1.0 Le=O.0T 
R^=1.50 Pr=1.0 Le=0.01 
R'=2.00 Pr=1.0 Le=0.01 
R*=5.00 Pr=1.0 Le=0.01 
NiI T = sqrt(Le=0.01)'Nu s 

100 



Figure 6. Nus versus NuT — 1 plot. Black dots represent the theoretical solutions. Each line 
stands for a different i?a* . Along each line the stability parameter decreases. The red lines denote 
convective sub-critial simulations, where Rax < Ra cr i t . 
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Discussion 

Trampedach: Your simulation of semi-convection looked 2-D. Is that correct? Will it 
be extended to 3-D? 

Muthsam: Yes, the present models are 2D. We will turn to 3D in the forseeable future. 

Gough: The stated purpose of the final solution you presented was to determine whether 
the interface between turbulent convection layer is purely diffusive, but you did not tell 
us whether that is the case. Please would you tell us whether the transport of heat and 
helium across the interface which you showed us, which maintains its integrity despite be- 
ing modulatory is purely diffusive or whether there is significant (small scale) convective 
transport? 

Muthsam: You obviously refer to the movie with the prc-assigncd step in the middle of 
the domain. The interface itself looks quite stable. Only on isolated spots (upstreams or 
downstreams) these streams can penetrate it somewhat. 

Pouquet: Do you see any advantage to use a hybrid (MPI/OpenMP) parallclization 
method? 

Muthsam: For really large runs you have to use MPI anyway for lack of shared memory. 
If you have, on one node, multiple cores you save messages when using them in the 
OpcnMP rather than in the MPI mode. 



